Liquid Limits: The Glass Transition and Liquid-Gas Spinodal Boundaries of 

Metastable Liquids 
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The liquid-gas spinodal and the glass transition define ultimate boundaries beyond which sub- 
stances cannot exist as (stable or metastable) liquids. The relation between these limits is analyzed 
via computer simulations of a model liquid. The results obtained indicate that the liquid - gas spin- 
odal and the glass transition lines intersect at a finite temperature, implying a glass - gas mechanical 
instability locus at low temperatures. The glass transition lines obtained by thermodynamic and 
dynamic criteria agree very well with each other. 
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Substances can exist in the liquid state beyond equi- 
librium phase boundaries in a metastable state, if nuclc- 
ation of the stable phase is avoided. Such metastable 
liquids are nevertheless bounded by ultimate limits in 
the form of the liquid-gas spinodal and the glass transi- 
tion locus. The liquid-gas spinodal is a limit of stability, 
beyond which the system is mechanically unstable. Al- 
though the nature of the glass transition, in particular 
the existence of an ideal glass transition, is still a matter 
of debate Q , the locus of glass transition temperatures 
defines a boundary beyond which substances transform 
to an amorphous solid state, and can no longer exist as 
liquids. A preliminary report of investigations concern- 
ing the relationship between these two limits of the liquid 
state are presented in this letter. 
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FIG. 1. (a) Isotherms from REMC simulations (points con- 
nected with lines), (b) Curves marked 'IS:' are pressure vs. 
density curves for inherent structures, with starting temper- 
atures T = 1.0,0.5. (c) Continuous curves labeled 'Eq. (3)' 
are from the empirical equation of state described in the text. 

The immediate motivation for the present study are 
observations made in Ref. B, where a threshold den- 



sity p* — 0.89 was identified for a model (mono) atomic 
liquid, across which the structure of typical local poten- 
tial energy minima or 'inherent structures' (IS) under- 
goes a qualitative change from spatially heterogeneous 
(below p*) to homogeneous structures (above p*). At 
the same density p* , the pressure P vs. density curve for 
the inherent structures goes through a minimum. Two 
speculations were made in [Q, motivated by these ob- 
servations: (1) The threshold density p* is (or closely 
approximates) the T — > limit of the liquid-gas spinodal 
locus p s (T), and (2) The threshold density p* forms the 
absolute lower density limit to glass formation. In the 
simplest scenario, the glass transition temperature ap- 
proaches zero as p — > p* from above. These speculations 
are tested here via computer simulations of a model liq- 
uid. The results obtained also prove valuable in assessing 
recent approaches to studying the glass transition. 

The model studied is a binary mixture of 204 type A 
and 52 type B particles, interacting via the Lennard- 
Jones (LJ) potential, with parameters zab/^aa = 1-5, 
^bb/^aa = 0.5, oabIoaa = 0.8, and o B bIoaa = 0.88, 
and ms/niA = 1- This system has been extensively 
studied as a model glass former [^|-^|. Monte Carlo 
(MC) simulations in the restricted ensemble (REMC) 
were performed (details below) for 8 temperatures for 
an average of 10 densities, for three different sets of con- 
straints in each case, to locate the liquid-gas spinodal 
locus. Run lengths ranged from 3 x 10 5 to 1.8 x 10 6 
MC cycles. Molecular Dynamics (MD) simulations (de- 
tails as in Q) were performed at 7 densities, p — 
1.08, 1.1, 1.125, 1.15, 1.2, 1.25, 1.35, over a wide range of 
temperatures (from T = 3.0 to 0.259 for p = 1.08). Run 
lengths ranged from 2 x 10 5 for T > 1.0 to 6.2 x 10 7 time 
steps (or 1.86 x 10 5 LJ time units, or 0.4 ps in Argon 
units). Simulations above T = 1.0 were typically done at 
constant temperature (NVT) and those below T — 1.0 at 
constant energy (NVE). NVT simulations were also per- 
formed for lower densities (p = 0.55 to 1.05, for T = 1.0 
to 0.325) to obtain isothermal compressibilities, kr- 

The liquid-gas spinodal is estimated via restricted en- 
semble simulations wherein the system is divided 
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into cells and fluctuations in density in each cell are re- 
stricted. Isotherms so obtained display a van der Waals- 
type loop, and permit estimation of the spinodal density 
from the location of their minima (See || for further de- 
tails) H . The resulting isotherms are displayed in Fig. 
1, along with the P vs. p curves for inherent struc- 
tures, obtained at two temperatures. The latter show 
that IS properties depend on the starting T, but the den- 
sity at the minimum is not significantly altered, and is 
p* = 1.08. Finite T isotherms do not show significant 
dependence on constraint strength. Spinodal densities 
p s (T) obtained from the location of minima along these 
isotherms, for each constraint strength, are shown (as 
spinodal temperatures T s (rho)) in Fig. 2. Also shown 
is p* where the IS pressure is minimum. To confirm the 
reliability of the REMC estimates, the spinodal is also 
estimated (a) by fitting isotherms (P vs. p) by (cubic) 
polynomials, and (b) by calculating kx directly in MD 
simulations, and locating the density where k^ 1 vanishes 
by polynomial extrapolation. The results are shown in 
Fig. 2, which shows that spinodal densities so obtained 
agree well with the REMC estimates. 
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FIG. 2. Liquid-gas spinodal obtained from (a) REMC sim- 
ulations, (b) kr, (c) polynomial fits to isotherms, and (d) the 
empirical free energy (T s Therm. Int.'). The same curve is 
also shown shifted in p by 0.08 ('T a (p - .08) Therm. Int.'). 
The glass transition locus obtained from (e) VFT fits to diffu- 
sivity data, and (d) extrapolation of configurational entropy 
to zero ('Tig Therm. Int.'). Also marked (*) is the density 
p* where inherent structure pressure is a minimum. 

An empirical free energy is constructed next, based on 
equilibrium data from MD simulations P,p|,p^|-p2[ , which 
is used both to obtain an independent estimate of the 
spinodal, and to obtain a thermodynamic estimate of the 
glass transition locus. The absolute free energy A(p, T) 
of the system at density p at a reference temperature 
T r = 3.0 is first defined in terms of the ideal gas con- 



tribution Aid (p, T) and the excess free energy A ex (p, T) 
obtained by integrating the pressure from simulations: 

A(p,T) = A id (p,T) + A ex (p,T), (1) 
(3A id (p,T) = N (3 InA + lnp-l), 



p r A ex (p,T r ) = p r A° ex (0,T r ) + N 
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Here, N is the number of particles, f3 = ksT, A is the de 
Broglie wavelength, and A° ex arises from the mixing en- 
tropy. A ex is fit to a fifth order polynomial is p Jl3| . A ex 
at a desired temperature may be evaluated by integrating 
the potential energy, E: 

f P > > 

(3A ex (p, (3) = l3A ex (p, p r )+ E(p, p )d(3 (2) 

As observed in HQ , the T dependence of E at fixed, high, 
density is well described by the form E(p,T) ~ T 3 / 5 , in 
agreement with predictions for dense liquids Jl4| . This 
form is not, however, accurate at low densities. In or- 
der to fit data well at low densities while retaining re- 
liability at high densities, E data for p = 0.55 to 1.35, 
and temperatures below T = 3.0 are fitted to the form 
E(p,T)/N = E (p) + E!(p)T E2 (P\ The parameters 
Eq,Ei,E2 are fitted polynomials in p p3[ . These fits, 
which represent the measured MD data very well, to- 
gether with the fit for A° x (0,T r = 3.0), define an empir- 
ical free energy from which the equation of state and the 
liquid-gas spinodal locus p s (T) are obtained via: 



P(P,T) = 



P 2 dA(p,T) 
N dp 



dP 
dp 



= 0. 



(3) 



The liquid-gas spinodal locus resulting from the em- 
pirical free energy is shown in Fig. 2 and is seen to oc- 
cur at lower densities than the REMC estimate, with a 
roughly constant shift in density (also shown in Fig. 2). 
Inspection of isotherms obtained from the empirical free 
energy (which show mean field behavior near the spin- 
odal) reveals that while they agree very well with simu- 
lation isotherms away from the spinodal, they show devi- 
ations very close to the (REMC) spinodal densities. Cor- 
respondingly fc^ 1 obtained in simulations drops faster to 
zero close to the spinodal, while at higher densities there 
is very good agreement. The estimate from the empirical 
free energy must thus be considered a lower limit to the 
location of the spinodal. 

Glass transition temperatures at the 7 studied densi- 
ties are obtained next, by fitting diffusion coefficients of 
A particles to the Vogel-Fulcher-Tammann-Hesse (VFT) 
form: 



D(p,T) = D (p)ex P 



A(p) 



T-T (p) 



(4) 
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Values of Tb(p) so estimated define lower limits to the 
laboratory glass transition T g . The diffusion coefficients 
and the corresponding VFT fits are shown in Fig. 3. The 
resulting glass transition temperatures To(p) are shown 
in Fig. 2. (Relaxation times yield very similar estimates. 
A qualitatively similar T g vs. p curve has recently been 
reported in |^5[.) 
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FIG. 3. Arrhenius plot of diffusion coefficients and the 
corresponding VFT fits for the 7 densities studied, (inset) 
Adam-Gibbs plot of diffusion coefficients against (TSc) -1 and 
straight line fits. 

A thermodynamic estimate of ideal glass transition 
temperatures is obtained next, by evaluating the con- 
figurational entropy of the liquid as a function of density 
and temperature. A number of studies in recent years 
P,|6|JfO|- |f^ , [l^ -^f| have focussed on evaluating of the con- 
figurational entropy [p2f . The approach here, based on 
analyzing inherent structures j23|, follows most closely 
those of Ref.s [5|j20|J^,p^| . The canonical partition func- 
tion is re- written as a sum over all local potential energy 
minima, which introduces a distribution function for the 
number of minima at a given energy: 

Q N (p,T)=A- 3N 1 ^ w - i Jdr N exp(-0$) (5) 
= V exp (-/3$ a ) A~ 3N [ dr N exp (-/3($ - 

= J d$ Q fi($ a ) exp(-(3($ a +Nf basin ($ a ,T))) 

= J d<S> Q exp(-j3($ a + Nf basin {<S> a ,T)-TS c ($ a ))) 

where $ is the total potential energy of the system, a in- 
dexes individual inherent structures, $ a is the potential 
energy at the minimum, f2($ Q ) is the number density 
of inherent structures with energy <!> Q , and the config- 
urational entropy S c = fcslnfi. The basin free energy 
fbasin(&a,T) is obtained by a restricted partition func- 



tion sum over a given inherent structure basin, V a . In 
the harmonic approximation, we have 

3N-3 

131 = 2 H h ] + 2iV^ lnAlE Pftherm + Pfmbl (6) 

i 

where are eigenvalues of the Hessian or curvature ma- 
trix at the minimum. f3f V i b is a slowly varying func- 
tion of temperature (the temperature dependence is ob- 
tained by averaging over 1 000, 100 inherent structures at 
T < 1., T > 1. respectively), and is fitted to the form 
Pfvibi.P,T) = fo(p) + fi(p)/T 2 which fits available data 
quite well. Polynomial fits to p-dependent parameters 
/o , /i are obtained |l3|] . The total entropy of the liquid 
S as well as the basin entropy Sbasin may be evaluated as 
a function of density and temperature from the total and 
basin free energies. The configurational entropy S c (p,T) 
and the ideal glass transition Tjg{p) are then given by, 

S c (p,T) = S(p,T) - S ba3in (p,T);S c {p,T IG (p)) = 0. (7) 

The ideal glass transition locus obtained is shown in Fig. 
2, and is seen to be in very good agreement with the es- 
timate based on diffusivity. Although such correlation is 
recorded for experimental data |24|] , it is noteworthy that 
To estimates here from dynamic data at quite high T pro- 
duce such agreement. This agreement offers mutual sup- 
port to the thermodynamic method followed here, and 
the feasibility of estimation based on dynamical measure- 
ments. T]c(p = 1.2) = 0.2976 is in very good agreement 
with the values 0.297 in || and 0.31 in Q. Diffusion co- 
efficients plotted against (TSc)^ 1 (inset of Fig. 3) show 
that the Adam-Gibbs expectation Q \nD — (T5 C ) _1 is 
extremely well satisfied at all densities [pd| , p^|j24[ . 

Data in Fig. 2 clearly indicate that the liquid-gas spin- 
odal and glass transition loci would intersect at a finite 
temperature T t - 0.16 (based on REMC data; - 0.12 
based on the empirical free energy). As this intersection 
happens at p ~ p* — 1.08, the data do not support the 
speculation that as T — > 0, p s {T) — > p* or pia{T) — > p* ■ 
On the other hand, p* does (within the uncertainty of the 
data) form the lower limiting density for glass formation. 
Experimental data (see e. g. ^5|) for the P-dependence 
of T g indicate that T g (P) typically intersects the zero 
pressure axis, implying glass - gas coexistence and finite 
T g at negative pressures, a possibility implicit in J2| and 
also noted in [^9). However, intersection of the liquid-gas 
spinodal and the glass transition locus at a finite temper- 
ature observed here further implies an ideal glass - gas 
mechanical instability below Ti . Indirect evidence for this 
exists in the form of experimental T g (P) loci displaying 
a negative slope ]2q ], but systematic evidence requires 
vitrification experiments at negative pressures p6[ . 

A proper study of such a mechanical instability is be- 
yond the scope of the present work, as it requires cal- 
culation of limits of stability accounting consistently for 
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the presence of the three relevant phases (liquid, gas, 
glass). As a preliminary attempt, the equation of state 
of the glass below Tjc(p) is obtained by noting that below 
Tig the system is trapped in the potential energy basin 
reached at Tig- Details of such a calculation will not be 
presented here; Instead, two consequences are presented 
in Fig. 4. The inset demonstrates the slope discontinuity 
of the pressure at fixed density at Tig for p = 1.2. The 
main diagram shows the isotherm at T = 0.05 for the 
glass phase, which displays a mechanical instability den- 
sity of p ~ 1., where the slope vanishes, or, kx diverges. 
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FIG. 4. Pressure vs. density for T = 0.05 in the ideal 
glass phase (solid line), which exhibits a minimum at p ~ 1.. 
Extension of the liquid's equation of state is shown for com- 
parison (dashed line). Inset shows pressure vs. temperature 
across the ideal glass transition for p — 1.2 (solid line). 

An interesting by-product of the present analysis is 
the estimation of configurational entropy as a function of 
both energy and density, and a range of temperature de- 
pendences of dynamical quantities displaying a variable 
degree of fragility |27| . These data permit an evaluation 
of the relationship between fragility and configurational 
entropy p£| , and will be presented elsewhere. 
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